Zero is not absence: censoring-based differential abundance analysis for microbiome data

Abstract Motivation Microbiome data analysis faces the challenge of sparsity, with many entries recorded as zeros. In differential abundance analysis, the presence of excessive zeros in data violates distributional assumptions and creates ties, leading to an increased risk of type I errors and reduced statistical power. Results We developed a novel normalization method, called censoring-based analysis of microbiome proportions (CAMP), for microbiome data by treating zeros as censored observations, transforming raw read counts into tie-free time-to-event-like data. This enables the use of survival analysis techniques, like the Cox proportional hazards model, for differential abundance analysis. Extensive simulations demonstrate that CAMP achieves proper type I error control and high power. Applying CAMP to a human gut microbiome dataset, we identify 60 new differentially abundant taxa across geographic locations, showcasing its usefulness. CAMP overcomes sparsity challenges, enabling improved statistical analysis and providing valuable insights into microbiome data in various contexts. Availability and implementation The R package is available at https://github.com/lapsumchan/CAMP.


S2 Sensitivity to choice of d
Building on the previous section's approach, we also utilized Simulation 1, with initial library size generated from Unif(25000, 40000) for both groups.However, we varied the value of d used in practice from d = 0.5 to d = 1, 2, 3, 4, and 5.
The choices of d = 0.5 and d = 1 are particularly noteworthy, as one reviewer highlighted d = 0.5 as an optimal choice from a measurement error perspective.
In Figure S1, we observe that the choice of d has a relatively minimal effect on type I error (Panel A).However, it significantly influences power (Panel B), with larger values of d leading to reduced power.This reduction in power may stem from the introduction of excessive noise when "imputing" non-zero entries with larger entries, potentially obscuring the distinction between nonnull and null taxa.Nevertheless, results for d = 0.5 and d = 1 are very similar, Censoring-based differential abundance analysis for microbiome data indicating no substantial practical difference between these values.From an interpretative standpoint, we recommend using d = 1, as it represents a more intuitive censoring cutoff for read counts, whereas d = 0.5 equivalent to half a read−lacks interpretability.

S3 Proportional hazards assumption
Both the log-rank test and Cox proportional hazards model (the case with covariates) are most powerful under the proportional hazards assumption.Following one reviewer's suggestion, we checked this assumption in the context of negative log-transformed relative abundance data.We used the same preprocessed human gut microbiome data from our real data analysis.For each taxa from the 3 pairwise comparison between countries, we tested the proportional hazards assumption using the Schonenfeld residuals.We applied the Benjamini-Hochberg procedure to the p-values of test and the proportion of Censoring-based differential abundance analysis for microbiome data 3 tests rejected at 5% level is presented in Table S1.As shown in Table S1, we failed to reject approximately 90% or more of the tests empirically, suggesting that there are minimal evidence against the proportional hazards assumption and potentially works fine as a working assumption.We hypothesize that there is a biological relevance as to why this assumption may hold: microbial communities can reach a state of relative equilibrium or balance such that the ratio of abundance or presence of different taxa remain relatively constant over abundance level.

S4 Performance metrics
We used type I error, power and false discover rate (FDR) to compare the performance of different DAA methods in simulations 1 to 3, where the number of true positive (TP) is the number of non-null taxa detected at the corresponding p-value threshold.Similarly, the number of true negatives (TN) is the number of null taxa not detected at the p-value threshold.Then, the number of false positives (FP) and the number of false negatives (FN) Censoring-based differential abundance analysis for microbiome data simply given by the number of non-null taxa minus TP, and the number of null taxa minus TN, respectively.For type I error, we used a 5% level; and we applied Benjamini Hochberg procedure to p-values at 5% level for power and FDR.Notice that FDR is not available for simulation 3 since the null and non-null taxa were generated separately.

S5 Real data analysis results
To complement our main simulations, which were based on data generated under assumed underlying models, we conducted a simulation using the same human gut microbiome data from the real data analysis.More information about this dataset can be found in Section 4 of Methods.After data preprocessing, we performed pairwise DAA by comparing the data from two countries at a time and repeating this process for all 3 2 = 3 country pairs.To assess the type I error rate, we permuted the country labels for each pair.
In this simulation (Figure S2), CAMP, LDM, and MetagenomeSeq were the only methods that consistently controlled the median type I error rate.
MetagenomeSeq, however, was slightly more conservative than the other two.
The remaining methods-corncob, DESeq2, and ZINQ-exhibited inflated median type I errors, specifically 20%, 12%, and 8%, respectively, across three pairwise comparisons.Once again, this showcases the advantages of employing a non-parametric approach in DAA.

S6 Sensitivity to library size
In response to a reviewer's comment that the negative log transformation depends on the library size, we explored the potential effects of differences in library sizes between groups.We modified Simulation 1 to generate the initial library size for group 1 from Unif(10000, 20000), while keeping the initial Censoring-based differential abundance analysis for microbiome data 9 library size for group 2 as Unif(25000, 40000), thus creating non-overlapping initial library sizes.We chose simulation 1, i.e., the Dirichlet-multinomial generation process without truncation, as it can capture the variability in microbial abundance due to biological factors more effectively.As observed in Figure S4, a systematic difference in library size can result in a slightly inflated type I error (Panel A) for CAMP, although the power (Panel B) remains relatively the same.
In practical settings, substantial and systematic differences in library size between groups are uncommon, barring technical variations like batch effects.
To address potential inflation in type I error due to library size differences, we propose a remedy.In cases where there is a statistically significant difference in library sizes between groups (e.g., from a t-test), one solution is to perform rarefaction before applying CAMP, thus, all the samples will have the exactly the same library size.

Fig. S1
Fig. S1 Median type I error and power for CAMP across 1000 replicates in simulation 1: (A) Type I error; (B) Power.

Fig. S2
Fig. S2 Type I error for 6 methods across 1000 replicates based on pairwise country comparison in simulation 4: (A) US vs Malawi; (B) Malawi vs Venezuela; (C): US vs Venezuela.

Fig. S3
Fig. S3 Differential abundance analysis results for gut microbiome dataset: (A) number of discoveries given by 4 methods (excluding corncob and DESeq2) in the Malawi vs US comparison; and (B) US vs Venezuela comparison.

Fig. S4
Fig. S4 Type I error and power for CAMP across 1000 replicates in a modified simulation 1: (A) Type I error; (B) Power.

Table S1
Proportion of tests that rejected the proportional hazards assumption after correcting for multiple comparison at 5% level in three pairwise comparisons (MA vs US; MA vs VE and US vs VE).

Table S2 :
List of 60 taxa (genus) uniquely identified by CAMP in the three pairwise comparisons (MA vs US; MA vs VE and US vs VE) with their corresponding CAMP p-values.p-values indicating statistical significance (at 5% false discovery rate) are bolded.